This document contains the code from “A comparison of statistical models used to characterize species-habitat associations with movement data”. This code includes the case study analysis on ringed seal movement in eastern Hudson Bay, Canada.
library(tidyverse)
library(ggplot2)
library(lubridate)
library(rnaturalearth)
library(rgdal)
library(terra)
library(amt)
library(momentuHMM)
library(viridis)
library(sf)
library(here)Next we load in the fish data. This is a subset of the data from Florko et al. (2021).
# load fish data
fish <- read.csv(here("data/preydiv.csv"))
head(fish)## lon lat preydiv
## 1 168329.9 -116109 0.513
## 2 168329.9 -166109 0.503
## 3 168329.9 -216109 0.512
## 4 168329.9 -266109 0.533
## 5 168329.9 -316109 0.557
## 6 168329.9 -366109 0.575
Visualize the fish data.
#prepare world data
## List of azimuthal equal area - HudsonBay
natearth <- ne_countries(scale = "medium",returnclass = "sp")
natearth <- natearth[which(natearth$continent!="Antarctica"),]
nat_trans <- spTransform(natearth,CRS("+proj=laea +lat_0=60 +lon_0=-85 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"))## Warning in spTransform(xSP, CRSobj, ...): NULL source CRS comment, falling back
## to PROJ string
## Warning in wkt(obj): CRS object has no comment
# plot fish map
fishmap <- ggplot() +
geom_tile(data = fish, aes(x = lon,y = lat,fill = preydiv)) +
scale_fill_viridis(option = "mako", limits = c(0.5, 0.75), name = "Prey\ndiversity") +
geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F)
fishmapNext we load in the seal movement data. This is a subset of the data from Florko et al. (2023).
seal <- read.csv(here("data/seal_track_m.csv"))
head(seal)## lon lat date id
## 1 357940.9 -368949.6 2012-10-29 14:00 116484_1
## 2 371594.3 -353092.9 2012-10-30 14:00 116484_1
## 3 388429.2 -361478.4 2012-10-31 14:00 116484_1
## 4 399326.2 -357288.1 2012-11-01 14:00 116484_1
## 5 411788.2 -349500.7 2012-11-02 14:00 116484_1
## 6 407708.1 -372564.6 2012-11-03 14:00 116484_1
# ensure the data is in the correct format
seal <- seal %>%
mutate(id = as.character(id),
date = as.Date(date)) Visualize the seal data on top of the fish data.
# plot
sealfishmap <- fishmap +
geom_point(data=seal, aes(x=lon, y=lat), alpha = 0.6, color = "#FCEEAE") +
geom_path(data=seal, aes(x=lon, y=lat), color = "#FCEEAE")
sealfishmapFit using amt (Signer et al. 2019)
# generate availability sample
set.seed(2023)
data_rsf <- seal %>%
make_track(lon, lat, date) %>%
random_points() # default is ten times as many available points as observed points
# plot used vs available locations
sealfishmap +
geom_point(data=data_rsf, aes(x=x_, y=y_, color = case_), alpha = 0.6) +
scale_color_manual(values = c("black", "#FCEEAE"),
label = c("Available", "Observed"), name = "Data type") # rasterize and extract prey diversity covariate
fish_raster <- terra::rast(fish, crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
data_rsf <- data_rsf %>%
extract_covariates(fish_raster)
# fit rsf
rsf1 <- data_rsf %>%
amt::fit_rsf(case_ ~ preydiv, model = TRUE)
# see parameter estimates
summary(rsf1)##
## Call:
## stats::glm(formula = formula, family = stats::binomial(link = "logit"),
## data = data, model = TRUE)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.5753 -0.4729 -0.4294 -0.3754 2.3363
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.118 1.400 -4.369 1.25e-05 ***
## preydiv 6.353 2.312 2.748 0.006 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 938.28 on 1539 degrees of freedom
## Residual deviance: 930.60 on 1538 degrees of freedom
## AIC: 934.6
##
## Number of Fisher Scoring iterations: 5
Fit using amt (Signer et al. 2019)
# generate availability sample
set.seed(2023)
data_ssf <- seal %>%
make_track(lon, lat, date) %>%
steps() %>%
random_steps() %>% # default is ten times as many available points as observed points
arrange(case_)
# plot used vs available locations
sealfishmap +
geom_point(data=data_ssf, aes(x=x2_, y=y2_, color = case_), alpha = 0.6) +
scale_color_manual(values = c("black", "#FCEEAE"),
label = c("Available", "Observed"), name = "Data type") # extract prey diversity covariate
data_ssf <- data_ssf %>%
extract_covariates(fish_raster, where = "end") #sample at end of step
# transform movement covariates
data_ssf <- data_ssf %>%
mutate(cos_ta = cos(ta_),
log_sl = log(sl_))
# fit ssfs
## ssf1: ssf without covariate affecting movement kernel
ssf1 <- data_ssf %>%
fit_clogit(case_ ~ preydiv + log_sl + cos_ta + strata(step_id_), model = TRUE)
## ssf2: ssf with covariate affecting movement kernel
ssf2 <- data_ssf %>%
fit_clogit(case_ ~ log_sl * preydiv + cos_ta + strata(step_id_), model = TRUE)
# see parameter estimates
summary(ssf1)## Call:
## coxph(formula = Surv(rep(1, 1518L), case_) ~ preydiv + log_sl +
## cos_ta + strata(step_id_), data = data, model = TRUE, method = "exact")
##
## n= 1516, number of events= 138
## (2 observations deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## preydiv 0.957799 2.605953 6.772711 0.141 0.888
## log_sl -0.008404 0.991631 0.083328 -0.101 0.920
## cos_ta 0.031707 1.032215 0.130643 0.243 0.808
##
## exp(coef) exp(-coef) lower .95 upper .95
## preydiv 2.6060 0.3837 4.477e-06 1.517e+06
## log_sl 0.9916 1.0084 8.422e-01 1.168e+00
## cos_ta 1.0322 0.9688 7.990e-01 1.333e+00
##
## Concordance= 0.494 (se = 0.029 )
## Likelihood ratio test= 0.09 on 3 df, p=1
## Wald test = 0.09 on 3 df, p=1
## Score (logrank) test = 0.09 on 3 df, p=1
summary(ssf2)## Call:
## coxph(formula = Surv(rep(1, 1518L), case_) ~ log_sl * preydiv +
## cos_ta + strata(step_id_), data = data, model = TRUE, method = "exact")
##
## n= 1516, number of events= 138
## (2 observations deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## log_sl 1.663e+00 5.273e+00 1.286e+00 1.293 0.196
## preydiv 2.791e+01 1.318e+12 2.177e+01 1.282 0.200
## cos_ta 2.338e-02 1.024e+00 1.311e-01 0.178 0.858
## log_sl:preydiv -2.744e+00 6.433e-02 2.101e+00 -1.306 0.192
##
## exp(coef) exp(-coef) lower .95 upper .95
## log_sl 5.273e+00 1.897e-01 4.241e-01 6.556e+01
## preydiv 1.318e+12 7.589e-13 3.919e-07 4.430e+30
## cos_ta 1.024e+00 9.769e-01 7.917e-01 1.324e+00
## log_sl:preydiv 6.433e-02 1.555e+01 1.046e-03 3.955e+00
##
## Concordance= 0.555 (se = 0.029 )
## Likelihood ratio test= 1.81 on 4 df, p=0.8
## Wald test = 1.81 on 4 df, p=0.8
## Score (logrank) test = 1.81 on 4 df, p=0.8
Fit using momentuHMM (McClintock and Michelot, 2018)
# prepare dataset
data_hmm <- seal %>%
make_track(lon, lat, date) %>%
extract_covariates(fish_raster) %>%
mutate(ID = 1,
x = x_,
y = y_,
date = t_) %>%
dplyr::select(ID, x, y, date, preydiv) %>%
as.data.frame()
# convert into momentuHMM format
data_hmm <- momentuHMM::prepData(data_hmm, coordNames = c("x", "y"),
type = "UTM",
covNames = c("preydiv"))
# define initial parameters
nbStates <- 3 # number of states
stepDist <- "gamma" # step distribution
angleDist <- "vm" # turning angle distribution
mu0 <- c(5000, 18000, 38000) # mean step length for each state
sigma0 <- c(5000, 10000, 10000) # sd step length for each state
stepPar0 <- c(mu0, sigma0)
kappa0 <- c(0.35, 0.55, 0.5) # turning angle for each state
formula = ~ preydiv # identify covariates
# fit basic 3-state hmm
set.seed(2023)
hmm1 <- momentuHMM::fitHMM(data=data_hmm, nbStates=nbStates,
dist=list(step=stepDist,angle=angleDist),
Par0=list(step=stepPar0,angle=kappa0))
# retrieve parameters to refine model
Par0_hmm1 <- momentuHMM::getPar0(hmm1, formula=formula)
# fit a refined hmm with parameters from hhm1
set.seed(2023)
hmm2 <- momentuHMM::fitHMM(data=data_hmm, nbStates=3,
dist=list(step=stepDist,angle=angleDist),
Par0=Par0_hmm1$Par,
delta0 = Par0_hmm1$delta,
beta0 = Par0_hmm1$beta,
formula=formula)
# plot decoded states
data_hmm$state <- as.factor(momentuHMM::viterbi(hmm2))
hmmstate_plot <- ggplot() +
scale_fill_viridis(option = "mako", limits = c(0.5, 0.75), name = "Prey\ndiversity") +
geom_path(data=data_hmm, aes(x=x, y=y, color = state, group =ID)) +
geom_point(data=data_hmm, aes(x=x, y=y, color = state, shape = state), size=2, alpha = 0.8) +
scale_color_manual(values = c("#99DDB6", "#539D9C", "#312C66"),
labels = c("Localized\nsearch", "Area-restricted\nsearch", "Travelling"),
name = "HMM state") +
scale_shape_manual(values = c(15, 16, 17),
labels=c("Localized\nsearch", "Area-restricted\nsearch", "Travelling"),
name = "HMM state") +
geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F)
hmmstate_plot#---------------------PLOT HISTOGRAMS
## prep data in lat/lon and refit model (so step length is in kilometers)
data_hmm <- seal %>%
make_track(lon, lat, date) %>%
extract_covariates(fish_raster) %>%
mutate(ID = 1,
x = x_,
y = y_,
date = t_) %>%
dplyr::select(ID, x, y, date, preydiv) %>%
as.data.frame() %>%
st_as_sf(coords = c("x", "y")) %>%
st_set_crs("+proj=laea +lat_0=60 +lon_0=-85 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0") %>%
st_transform("+proj=longlat +datum=WGS84") %>%
mutate(long = unlist(map(geometry,1)),
lat = unlist(map(geometry,2))) %>%
st_drop_geometry()
# do momentuhmm which calculates step length in KM
data_hmm <- momentuHMM::prepData(data_hmm, coordNames = c("long", "lat"), type = "LL", covNames = c("preydiv"))
# define parameters
nbStates <- 3
stepDist <- "gamma" # step distribution
angleDist <- "vm" # turning angle distribution
mu0 <- c(5, 12, 38)
sigma0 <- c(3, 5, 8)
stepPar0 <- c(mu0, sigma0)
kappa0 <- c(0.35, 0.55, 0.5)
# fit HMM (with step in km)
set.seed(2023)
hmm1_km <- momentuHMM::fitHMM(data=data_hmm, nbStates=nbStates,
dist=list(step=stepDist,angle=angleDist),
Par0=list(step=stepPar0,angle=kappa0))
# get parameters
Par0_hmm1_km <- momentuHMM::getPar0(hmm1_km, formula=formula)
# fit HMM with parameters from hmm1_km
set.seed(2023)
hmm2_km <- momentuHMM::fitHMM(data=data_hmm, nbStates=3,
dist=list(step=stepDist,angle=angleDist),
Par0=Par0_hmm1_km$Par,
delta0 = Par0_hmm1_km$delta,
beta0 = Par0_hmm1_km$beta,
formula=formula)
# add the state estimate from the HMM
data_hmm$state <- as.factor(momentuHMM::viterbi(hmm2_km))
##-----STEP DIST histogram
#calculate frequencies of states
v <- momentuHMM::viterbi(hmm2_km)
stateFreq <- table(v) / length(v)
#plot colours
colours.states <- c("#99DDB6", "#539D9C", "#312C66")
#generate sequence for x axis of density functions
x <- seq(0, 50, length=1000)
#get converged mean and sd for each state
meanARS <- hmm2_km$mle$step[1,1]
sdARS <- hmm2_km$mle$step[2,1]
meanCR <- hmm2_km$mle$step[1,2]
sdCR <- hmm2_km$mle$step[2,2]
meanTR <- hmm2_km$mle$step[1,3]
sdTR <- hmm2_km$mle$step[2,3]
#calculate shape and scale of the gamma distributions from mean and sd
sh <- function(mean, sd) { return(mean^2 / sd^2)}
sc <- function(mean, sd) { return(sd^2 / mean)}
#get density functions of the distributions
y_ARS <- dgamma(x, shape=sh(meanARS,sdARS), scale=sc(meanARS,sdARS)) * stateFreq[[1]]
y_CR <- dgamma(x, shape=sh(meanCR,sdCR), scale=sc(meanCR,sdCR)) * stateFreq[[2]]
y_TR <- dgamma(x, shape=sh(meanTR,sdTR), scale=sc(meanTR,sdTR)) * stateFreq[[3]]
#combine densities in a single dataframe for more convenient plotting
df.y_ARS <- data.frame(dens=y_ARS, State="Foraging", x=x)
df.y_CR <- data.frame(dens=y_CR, State="ARS", x=x)
df.y_TR <- data.frame(dens=y_TR, State="Travelling", x=x)
statedis <- rbind(df.y_ARS, df.y_CR, df.y_TR)
#plot distributions
hmmstepdist_plot <- ggplot() +
geom_line(data=statedis,aes(x=x,y=dens,colour=State,linetype=State), size=1.2) +
scale_colour_manual(values=c(colours.states,"#000000"),
breaks = c('Foraging', 'ARS', 'Travelling'),
labels=c("Localized\nsearch", "Area-restricted\nsearch", "Travelling")) +
scale_linetype_manual(values=c("solid","solid", "solid"),
breaks = c('Foraging', 'ARS', 'Travelling'),
labels=c("Localized\nsearch", "Area-restricted\nsearch", "Travelling")) +
ylab("Density") +
xlab("Step length (kms)") +
theme_minimal()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## Warning: Please use `linewidth` instead.
hmmstepdist_plot##-----ANGLE DIST
# #generate sequence for x axis of density functions
x <- seq(-pi, pi,length=1000)
#get converged mean and concentration for each state
meanARS <- hmm2_km$mle$angle[1,1]
sdARS <- hmm2_km$mle$angle[2,1]
meanCR <- hmm2_km$mle$angle[1,2]
sdCR <- hmm2_km$mle$angle[2,2]
meanTR <- hmm2_km$mle$angle[1,3]
sdTR <- hmm2_km$mle$angle[2,3]
#get density functions of the distributions
y_ARS <- CircStats::dvm(x, mu=meanARS, kappa=sdARS) * stateFreq[[1]]
y_CR <- CircStats::dvm(x, mu=meanCR, kappa=sdCR) * stateFreq[[2]]
y_TR <- CircStats::dvm(x, mu=meanTR, kappa=sdTR) * stateFreq[[3]]
#combine densities in a single dataframe for more convenient plotting
df.y_ARS <- data.frame(dens=y_ARS, State="Foraging",x=x)
df.y_CR <- data.frame(dens=y_CR, State="ARS",x=x)
df.y_TR <- data.frame(dens=y_TR, State="Travelling",x=x)
cmb <- rbind(df.y_TR, df.y_CR, df.y_ARS)
#plot distributions
hmmangledist_plot <- ggplot() +
geom_line(data=cmb,aes(x=x,y=dens,colour=State), size = 1.2) +
scale_colour_manual(values=c(colours.states), breaks = c('Foraging', 'ARS', 'Travelling'), labels=c("Localized\nsearch", "Area-restricted\nsearch", "Travelling")) +
scale_x_continuous(limits=c(-pi,pi))+
ylab("Density") +
xlab("Turn angle (radians) ") +
theme_minimal()
hmmangledist_plotPrep the data
# prep the fish data
newfish <- fish_raster %>%
terra::as.data.frame(xy = TRUE) %>%
filter(x > 100000 & x < 600000 & y > -550000 & y < 0)modcoef <- summary(rsf1)$coef
x <- exp(modcoef[1] + (modcoef[2] * newfish$preydiv))
x <- x / (1 + x)
# range fn
range01 <- function(x){(x-min(x))/(max(x)-min(x))}
newfish$rsf_prediction <-range01(x)
# plot
map_RSF <- ggplot() +
geom_tile(data = newfish, aes(x = x,y = y, fill = rsf_prediction)) +
scale_fill_viridis(option = "mako", name = "RSF prediction", limits = c(0, 0.3)) +
geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F) ## Regions defined for each Polygons
map_RSF# generate availability sample
set.seed(2023)
data_ssf <- seal %>%
mutate(date = as.POSIXct(date)) %>%
make_track(lon, lat, date) %>%
steps() %>%
random_steps() %>%
arrange(case_) %>%
amt::extract_covariates(fish_raster, where = "both") %>%
na.omit()
### -------------- SSF1
# fit SSF1 model
m1 <- data_ssf |>
fit_clogit(case_ ~ preydiv_end + cos(ta_) + log(sl_) +
# add terms for home range
x2_ + y2_ + I(x2_^2 + y2_^2) +
strata(step_id_))
# Start with simulations
# First generate a redistribution kernel
set.seed(2023)
start <- make_start((seal %>%
mutate(date = as.POSIXct(date)) %>%
make_track(lon, lat, date))
# random sample along the track as the start
[sample(1:nrow(seal), 1), ],
dt_ = hours(2))
# generate redistribution kernel
k1 <- redistribution_kernel(m1, map = fish_raster, start = start,
stochastic = TRUE,
tolerance.outside = 0.1,
n.control = 1e3)
# Now simulate a path of length 30 for 50 times. This will lead a transient UD. For a steady state UD, either increase `n` or choose many more random starting points.
n_steps = 1e3
n_steps1 = n_steps + 1
set.seed(2023)
p1 <- amt::simulate_path(x = k1, n = n_steps, start = start, verbose = TRUE)
# plot simulated track
ssf_track_1 <- fishmap +
geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
geom_point(data = p1, aes(x = x_,y = y_), alpha = 0.61) +
geom_path(data = p1, aes(x = x_,y = y_)) ## Regions defined for each Polygons
# estimate SSUD
uds_ssf1 <- tibble(rep = 1:n_steps1,
x_ = p1$x_, y_ = p1$y_,
t_ = p1$t_, dt = p1$dt) |>
filter(!is.na(x_)) |>
make_track(x_, y_) |>
hr_kde(trast = fish_raster, which_min = "global") %>%
hr_ud() %>%
terra::as.data.frame(xy = TRUE)
# plot SSUD
map_SSF1 <- ggplot() +
geom_tile(data = uds_ssf1, aes(x = x,y = y, fill = preydiv)) +
scale_fill_viridis(option = "mako", name = "SSF1 prediction") +
geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F) ## Regions defined for each Polygons
map_SSF1### -------------- SSF2
# fit SSF2 model
m2 <- data_ssf |>
fit_clogit(case_ ~ preydiv_end + preydiv_end:cos(ta_)+
preydiv_end:log(sl_) +
# add terms for home range
x2_ + y2_ + I(x2_^2 + y2_^2) +
strata(step_id_))
# Start with simulations
# First generate a redistribution kernel
set.seed(2023)
start <- make_start((seal %>%
mutate(date = as.POSIXct(date)) %>%
make_track(lon, lat, date))
# first location of the track as the start
[sample(1:nrow(seal), 1), ],
dt_ = hours(2))
# generate redistribution kernel
k2 <- redistribution_kernel(m2, map = fish_raster, start = start,
stochastic = TRUE,
tolerance.outside = 0.3,
n.control = 1e3)
# Now simulate a path of length 30 for 50 times. This will lead a transient UD. For a steady state UD, either increase `n` or choose many more random starting points.
n_steps = 1e3
n_steps1 = n_steps + 1
set.seed(2023)
p2 <- amt::simulate_path(x = k2, n = n_steps, start = start)
# plot simulated track
ssf_track_2 <- fishmap +
geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
geom_point(data = p2, aes(x = x_,y = y_), alpha = 0.61) +
geom_path(data = p2, aes(x = x_,y = y_)) ## Regions defined for each Polygons
# estimate SSUD
uds_ssf2 <- tibble(rep = 1:n_steps1,
x_ = p1$x_, y_ = p1$y_,
t_ = p1$t_, dt = p1$dt) |>
filter(!is.na(x_)) |>
make_track(x_, y_) |>
hr_kde(trast = fish_raster, which_min = "local") %>%
hr_ud() %>%
terra::as.data.frame(xy = TRUE)
# plot SSUD
map_SSF2 <- ggplot() +
geom_tile(data = uds_ssf2, aes(x = x,y = y, fill = preydiv)) +
scale_fill_viridis(option = "mako", name = "SSF2 prediction") +
geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F) ## Regions defined for each Polygons
map_SSF2library(cowplot)##
## Attaching package: 'cowplot'
##
## The following object is masked from 'package:lubridate':
##
## stamp
plot_grid(map_SSF1, map_SSF2, ssf_track_1, ssf_track_2, nrow =2)x <- as.data.frame(momentuHMM::stationary(hmm2, data.frame(preydiv = newfish$preydiv)))
newfish$hmm_state1 <- x$state.1
newfish$hmm_state2 <- x$state.2
newfish$hmm_state3 <- x$state.3
newfish_long <- newfish %>%
tidyr::pivot_longer(cols = hmm_state1:hmm_state3,
names_to = "model", values_to = "prediction") %>%
mutate(dplyr::across(model, factor, levels=
c("hmm_state1", "hmm_state2", "hmm_state3")))
map_HMM <- ggplot() +
geom_tile(data = newfish_long, aes(x = x, y = y, fill = prediction)) +
scale_fill_viridis(option = "mako", limits = c(0,1)) +
labs(fill = 'HMM predicted\nprobability') +
geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), colour = "white", fill = "grey90", size = 0.3) +
coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F) +
facet_wrap(~model, labeller = as_labeller(c('hmm_state1' = "Localized search",
'hmm_state2' = "Area-restricted search",
'hmm_state3' = "Travel")))## Regions defined for each Polygons
map_HMMbase <- newfish %>%
# (note, it doesn't matter what values you pick for sl_ and ta_)
mutate(log_sl = log(45),
cos_ta = cos(1))
x1 <- base
x2 <- base %>%
mutate(preydiv = mean(base$preydiv))
# So how do we iterate this? We could use a for loop, we could use lapply, etc.
# I will use lapply here, with a custom function inside.
log_rss_list <- lapply(1:nrow(x1), function(i) {
# Report status
cat(i, "of", nrow(x1), "\n")
# Calculate log-RSS for that row
xx <- log_rss(rsf1, x1[i,], x2[i,], ci = "se")
# Return the element $df
return(xx$df)
})## 1 of 72
## 2 of 72
## 3 of 72
## 4 of 72
## 5 of 72
## 6 of 72
## 7 of 72
## 8 of 72
## 9 of 72
## 10 of 72
## 11 of 72
## 12 of 72
## 13 of 72
## 14 of 72
## 15 of 72
## 16 of 72
## 17 of 72
## 18 of 72
## 19 of 72
## 20 of 72
## 21 of 72
## 22 of 72
## 23 of 72
## 24 of 72
## 25 of 72
## 26 of 72
## 27 of 72
## 28 of 72
## 29 of 72
## 30 of 72
## 31 of 72
## 32 of 72
## 33 of 72
## 34 of 72
## 35 of 72
## 36 of 72
## 37 of 72
## 38 of 72
## 39 of 72
## 40 of 72
## 41 of 72
## 42 of 72
## 43 of 72
## 44 of 72
## 45 of 72
## 46 of 72
## 47 of 72
## 48 of 72
## 49 of 72
## 50 of 72
## 51 of 72
## 52 of 72
## 53 of 72
## 54 of 72
## 55 of 72
## 56 of 72
## 57 of 72
## 58 of 72
## 59 of 72
## 60 of 72
## 61 of 72
## 62 of 72
## 63 of 72
## 64 of 72
## 65 of 72
## 66 of 72
## 67 of 72
## 68 of 72
## 69 of 72
## 70 of 72
## 71 of 72
## 72 of 72
# We now have a list with an element for each row. We can combine those rows
# into a single data.frame with 'dplyr::bind_rows()'
res <- dplyr::bind_rows(log_rss_list)
# PLOT
rsf_line <- ggplot(res, aes(x = preydiv_x1, y = (log_rss))) +
geom_line(size = 1, color = "#F8D59F",linetype = 2) +
geom_ribbon(aes(ymin=lwr, ymax=upr, x=preydiv_x1), alpha = 0.4, fill = "#F8D59F") +
xlab("Prey diversity") +
ylab("log-RSS") +
theme_minimal()## Prediction for SSF1
x1 <- base
x2 <- base %>%
mutate(preydiv = mean(base$preydiv))
# So how do we iterate this? We could use a for loop, we could use lapply, etc.
# I will use lapply here, with a custom function inside.
log_rss_list <- lapply(1:nrow(x1), function(i) {
# Calculate log-RSS for that row
xx <- log_rss(ssf1, x1[i,], x2[i,], ci = "se")
# Return the element $df
return(xx$df)
})
# We now have a list with an element for each row. We can combine those rows
# into a single data.frame with 'dplyr::bind_rows()'
res <- dplyr::bind_rows(log_rss_list)
# PLOT
ssf_line_1 <- ggplot() +
geom_line(data = res, aes(x = preydiv_x1, y = log_rss), size = 1, linetype = 3) +
geom_ribbon(data = res, aes(ymin=lwr, ymax=upr, x=preydiv_x1), alpha = 0.3) +
xlab("Prey diversity") +
ylab("log-RSS") +
theme_minimal()
## Prediction for SSF2
nums = seal %>%
make_track(lon, lat, date) %>%
steps %>%
summarize(quants = quantile(sl_, c(0.25, 0.5, 0.75))) %>%
pull()
x1 <- base
results <- lapply(nums, function(i) {
x1$log_sl <- i
x2 <- base %>%
mutate(preydiv = mean(base$preydiv))
# So how do we iterate this? We could use a for loop, we could use lapply, etc.
# I will use lapply here, with a custom function inside.
log_rss_list <- lapply(1:nrow(x1), function(i) {
# Calculate log-RSS for that row
xx <- log_rss(ssf2, x1[i,], x2[i,], ci = "se")
# Return the element $df
return(xx$df)
})
# We now have a list with an element for each row. We can combine those rows
# into a single data.frame with 'dplyr::bind_rows()'
res <- dplyr::bind_rows(log_rss_list)
} )
results <- dplyr::bind_rows(results) %>%
mutate(log_sl_x1 = as.factor(round(log_sl_x1, 1)))
# PLOT
ggplot(results, aes(x = preydiv_x1, y = (log_rss))) +
geom_line(size = 1, aes(color = log_sl_x1, group = log_sl_x1)) +
geom_ribbon(aes(ymin=lwr, ymax=upr, x=preydiv_x1, fill = log_sl_x1, group = log_sl_x1), alpha = 0.3) +
xlab("Prey diversity") +
ylab("log-RSS") +
theme_minimal()# Plot
ssf_line_1 +
geom_line(data = results, size = 1, aes(x = preydiv_x1, y = (log_rss),
color = log_sl_x1, group = log_sl_x1)) +
geom_ribbon(data = results, aes(ymin=lwr, ymax=upr, x=preydiv_x1,
fill = log_sl_x1, group = log_sl_x1), alpha = 0.3)# grab stationary probabilities
ps <- momentuHMM::plotStationary(hmm2, plotCI= TRUE, return = TRUE)state1 <- ps$preydiv$'state 1' %>% mutate(state = 1)
state2 <- ps$preydiv$'state 2' %>% mutate(state = 2)
state3 <- ps$preydiv$'state 3' %>% mutate(state = 3)
# bind to one data frame
pdat <- rbind(state1, state2, state3) %>%
mutate(state = as.character(state)) # plot
hmmprobs_plot <- ggplot() +
geom_line(data = pdat, aes(x = cov, y = est, color = state)) +
geom_ribbon(data = pdat, aes(x=cov, y=est, ymax=est+se, ymin=est-se, fill = state),
alpha = 0.4, show.legend = TRUE) +
ylab("Stationary state probabilties") +
xlab("Prey diversity") +
scale_color_manual(values = c("#99DDB6", "#539D9C", "#312C66"), name = "HMM state",
labels=c("Localized\nsearch", "Area-restricted\nsearch", "Travel")) +
scale_fill_manual(values = c("#99DDB6", "#539D9C", "#312C66"), name = "HMM state",
labels=c("Localized\nsearch", "Area-restricted\nsearch", "Travel")) +
theme_minimal()
hmmprobs_plot